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Abstract 

We propose a new geometric method of IR factorization in sector de- 
composition. The problem is converted into a set of problems in convex 
geometry. The latter problems are solved using algorithms in combinato- 
rial geometry. This method provides a deterministic algorithm and never 
falls into an infinite loop. The number of resulting sectors depends on 
the algorithm of triangulation. Our test implementation shows smaller 
number of sectors comparing with other existing methods with iterations. 

1 Introduction 

Analysis of recent high statistical experiments in high energy physics requires 
accurate theoretical predictions. Higher order corrections in perturbation theory 
become indispensable for the calculation of many physical processes. Hovifever, 
calculation of higher loop integration is not an easy problem because of its high 
dimensionality and complicated dependence on many kinematic and mass pa- 
rameters. As the required final results of calculations are sets of numerical val- 
ues, various numerical methods of loop integrations have been proposed. Since 
numerical method can only be applied to the quantities without divergences, 
one has to avoid singularities in the integrand by regularization or subtraction 
of divergent part before applying numerical method. 

In the analysis of loop corrections in perturbative QCD, D-dimensional reg- 
ularization of infrared (IR) divergences is commonly used, where these diver- 
gences are separated in the form of poles in terms of e = (4 — D)/2. The 
method of sector decomposition [2 [H [31 [H El [5] provides a systematic proce- 
dure to extract these poles, and their coefficients are free from IR divergences. 
This method consists of several steps, whose details will be discussed in the next 
section. Here we mention each of these steps briefly: 

1. Start from Feynman parameter representation of a loop integration. 

2. Primary sector decomposition: 
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Divide the domain of integration into sectors, eliminating the (5-distribution, 
such that each sector includes one corner of the original integration do- 
main. 

3. Factorization step: 

Decompose sectors repeatedly such that IR diverging part of the integrand 
is factored out as a product of powers of the integration variables in each 
sector. 

4. Divergence separation step: 

When the integrand in each sector has a divergent factor of an integration 
variable, one can extract the singular part expanding into Laurent series. 
Poles of e are separated with the analytic integration of the singular part. 
After all poles are extracted by repeating this procedure for all divergent 
factors, the finite part is expanded into Taylor series in terms of e. 

5. Integration step: 

By the previous step, the coefficients of poles and the finite part are ex- 
pressed by integrals of functions without IR divergences. If these functions 
do not suffer from other singularities, these integrations are calculated nu- 
merically. Otherwise further regulations of other singularities are neces- 
sary [ZllHlinillD] 

The factorization step of the usual sector decomposition approach decom- 
poses the integration domain into sectors according to a sector decomposition 
strategy. A strategy provides a rule to choose a subset of integration variables 
and to determine how chosen variables are transformed and how the domain is 
divided. After applying these decompositions iteratively, one will obtain inte- 
grands whose IR divergent part is correctly factored out in each resulting sector. 
However, it is found that some strategies may lead to an infinite loop of iteration. 
This problem was first solved by Ref. [HI [12] with presenting three strategies 
guaranteed to terminate. Their strategies are based on the analysis of the poly- 
nomials appearing in the integrand. The strategy given in Ref. [SI [ID- also 
solve this problem, which is based on the analysis of polytopes corresponding 
to the polynomials. 

The number of generated sectors depends heavily on the strategies. The 
efficiency of sector decomposition method depends on the resulting number of 
decomposed sectors, since the numerical integration is the most time consuming 
part and will be repeated for different values of kinematic or mass parameters. 
Therefore such a method of the factorization step is favored that gives as small 
a number of generated sectors as possible. 

In this paper, we propose a new method of the factorization step employing 
geometric interpretation of the problem. This is a deterministic method without 
iterated decomposition and hence never falls into an infinite loop. In section 2, 
the problem is rephrased in clarifying the condition of the expected output of 
the factorization step. The problem is converted to a set of problems in convex 
geometry in section 3. Latter problems are solved in section 4 with algorithms 
developed in combinatorial geometry. Then it is shown that the number of the 
decomposed sectors depends on the triangulation algorithm of n-dimensional 
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polytopes. In section 5, we show the resuhs of the test hxiplementation of our 
method. The last section is devoted to the conclusion. 

Throughout this article, we use the following 'multi-index notation'. For m- 
tuple of variables x — {xi, • • • , Xm) and numbers a — (ai, • • • , am), the product 
of a;"^ for j = 1 , • • • , to is expressed by 

X := x^^ ■ ■ ■ x^ , (1) 

and 

d™x := dxi ■ ■ ■ dxm- (2) 
We also use special vector I defined byl:=(l,---,l). 

2 Sector decomposition 

Let us start with Eq. (11) of Ref. [S] for the loop integration expressed with 
Feynman parameters x = {xi, ■ • ■ , xn)- 

where 

• i is the number of loops, 

• _D = 4 — 2e is the dimension of the space-time, 

• is the dimension of the integration, 

• '^j — li"'' is the power of the propagator corresponding to the 
Feynman parameter x^ , 

• U\s& homogeneous polynomial of {xj} of degree L, and all the coefficients 
of the monomials of lA are equal to 1. 

• J-" is a homogeneous polynomial of {xj} of degree L+1, and the coefficients 
of the monomials of J- consist of kinematic and mass parameters. 

Polynomial may approach to when the vector x of Feynman param- 
eters approaches to a corner of the integration domain constrained by the S- 
distribution [xj except for one parameter). In such cases, the integration 
may diverge and cause IR singularities, which are dimensionally regulated. 

The primary sector decomposition step is to separate these singularities with 
dividing the integration domain into N sub-domains such that each sector in- 
cludes only one corner. Sector I is defined by: 

{{xi,X2, ■ ■ ■ ,Xn) I Xj < Xl, Vj ^ I}. (4) 
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In this sector, the integration variables are changed from {xi, • • • , xn} to {xi,ti, • • • , ijv-i} 
in the following way: 

Xj = xit'j, (5) 

{tj for j < I 
1 for j=l . (6) 
tj_i for j > I 

Integrating over xi in Eq. ([3]), one obtain Eq. (14) of Ref. [H]: 

nf=ir(..) 

Gi ^ [ d''-hr'~^u^{t):F^{t), (8) 

Jo 



where 



and 



Ti{t) = Tit'), (3 ^ ^{N,-LD/2), (9) 
Ui{t) = U{t'), ^ = N,-{L+l)D/2, (10) 



i^j for j < ^ 
^■j+i for j > I 



(11) 



After the primary sector decomposition, the domain of integration becomes (N— 
l)-dimensional unit cube and the corner of the original domain of integration is 
mapped to the origin tj = (Vj), where Ui and J-i may vanish. This expression 
is the starting point of our discussion. 

The purpose of the factorization step is factoring out the possible singular- 
ities at tj — from Ui and Ti. In the usual sector decomposition approach, 
this factorization is obtained by dividing the domain of integration repeatedly 
and applying appropriate transformation of integration variables. With an ap- 
propriate strategy, one can finally find proper sub-domains and new variables 
suitable for the separation of IR divergences. To be more specific: 

1. Integration domain is divided into m sectors {Di, D2, • • • , Dm}- 

2. In each sector Da, new variables z = (zi, • • ■ ,zjv-i) are introduced and 
tj = tj{z) is expressed as a monomial of z, and Jacobian Ja{z) is a mono- 
mial of z with constant coefficient. Let monomial z"^" be a product of 
and Jacobian: 

= n^f°^' •= t{-^T'~^Ja{z) - \[tk{zr'''-^Ja{z). (12) 



3. Integration domain of z is the (A^ — l)-dimensional unit cube (0 < Zj < 1). 

4. Polynomials Ui and J-i are expressed in the following form in each sector 

Da-. 

Ul = Caz"^ {1 + Ha{z)) , (13) 

= cy^ii+H'M}. (14) 
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where ba and b'^ arc {N — l)-tuplcs of non-negative integers {ba,b'^ £ 
Z^o~^fl, Ha{z) and H'^{z) are polynomials of z such that Ha{Q) = and 
Hl{Q) = 0. 

The output of the factorization step is shown by the following form of the 
expression: 



(15) 



The divergence separation step proceeds in the following way. Factor 2;Ca+baf'+''£i7 



is expressed by z^' ' {rij^mj £ Z). If < — 1, the integration of Zj di- 
verges when e ^ 0. In this case, the integrand is split into IR divergent and 
finite parts. Let us consider the case of Uj = —I for simplicity: 



dz.z-''-"''' il + Haiz)ril + Kiz)f 

' {1 + Ha{z)U^=Oy {I + Ha'{z)U^=o) 
[il+Ha{z)y{l+Ha>iz)f 
[l + Ha{z)U^=Oy (l+Ha'{z)L=0) 



dZj Zj 



dZj Zj 



(16) 



The first term is integrated to 
1 



{I + Ha{z)U^^Oy {I + Ha'iz^^^oy , 



(17) 



and the integrand of the second term is finite for Zj 0. When rij < —1, 
it is necessary to expand (1 -I- Ha{z)y (1 -|- Ha'{z)Y around Zj = up to the 
necessary order for the separation of the IR divergent parts. Repeating this 
procedure for all divergent integrals, we can extract all IR divergent parts as 
poles in terms of e. Then we can freely expand integrands of remaining integrals 
with respect to e, and obtain a Laurent series of the whole integral in terms of e, 
whose coefficients are expressed by IR finite integrals. We can calculate these IR 
finite integrals numerically if there are no other singularities in the integrands. 

Let us go back to the factorization step. In order to see the problem in a 
different way, let us apply the expression to a simple polynomial case of 
three variables ti,t2,t^'. 



When the integration variables are changed to zi, Z2, z^ by 



h = 2:12:223, 



^2 = 21, 



^3 = 22, 



(18) 



(19) 



^In this paper, the set of non- negative integers and the set of non-negative real numbers 
are denoted by Z>o,R>o, respectively. 
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we obtain 



Tl{t{z)) = 2:iZ2(-S23 - S4Z2Z3 - Sl2Zz). (20) 

This expression implies that around the origin of the domain {(^1(2:), i2(2;), t3(0)) | < 
Zi < 1}, the term 523^2^3 = S2'iZiZ2 dominates to the others and the ratios of 
other terms to the dominant one become monomials of the new variables zi. 
Similarly, changing variables in Eq. (|18p as 

ii=2i, t2 = ziZ2, h = Z3, (21) 

gives us 

Ti{t{z)) = zi{-si2 - S2ZZ2ZZ - 54^3), (22) 

and this exhibits a case that 512^1 = Si2-Zi is dominant at the origin. Note that 
54^1^3 cannot be a dominant term at the origin because always |si2ii| > 154^1^3! 
when all tj approach to 0. 

Above example shows that the factorization step is equivalent to solve the 
following set of problems: 

1. to determine which terms can be dominant, 

2. to determine the sub-domain where each term becomes dominant, and 

3. to find a new parameterization of the variables such that other terms are 
expressed as monomials multiplied by the dominant term. 

In the following sections, we will show that these problems can be solved by 
a new method without iterations. 



3 Geometric method 

Let us change the integration variables from t = {ii, • • ■ , ijv-i} to y = {t/i, • ■ • , y^-i) 
in Eq. ^ as the following: 

h^e'^^ or y, = -logi, (j = 1, • • • , iV - 1). (23) 

Jacobian is ^3 and the integration domain is changed from < < 1 to 
< Vj < 00. We obtain 

/"OC 

Gi = / d''-^ye~^''''y^U;'{e-y)T^{e-y), (24) 
Jo 

where {v' ,y) is the inner product defined in {N — l)-dimensional Euclidean 
space. 

A monomial = t^^ t^2 ' ' ' ^^n-i characterized by an integer vector h — 
{&i, • • • , fejv-i} {bj G Z>q"^). Since a polynomial P is a sum of monomials, the 
polynomial corresponds to a finite set of such vectors C Z>(7^- In terms of 
the new variables {y^}, P is expressed by: 

P{t) = J2 ^bt'' = J2 ate-'^'^'y^ (a& 7^ 0, Vfe G Z^). (25) 

bezp bezp 
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Let us consider the limit A 



for the variables yj = Xu 



-oo (Vj) 



with a constant vector u G R"q ^\{0}I1 The dominant term of P(t) in this limit 
is the monomial a^t'' whose value of {b, u) takes the minimum value of (c, u) 
among all c G Z^: 



{b,u) = min{(c, u) I c G Z"^}. 



(26) 



In other words, a monomial for a vector b is dominant in the domain of y 
space, which is defined by: 



{z/gK^o-1 I (6,y) < {c,y), VceZ^}. 



(27) 



Evidently the union of all covers whole integration domain. The polynomial 
P of Eq. dm is rewritten to: 



my)) - E ^(2^ ^ o ^"'"'"^ 



at 



cezp\{b} 



-{c-b.y) 



(28) 



Since {c — b,y) > for y G A^, the behavior of P in the limit of || y \\^ oo is 
determined by the term abe~^'''y\ If a monomial corresponding to b cannot be 
a dominant term, A^ is the empty set. Including this case, when the dimension 
of A^ is less than — 1, the {N — l)-dimcnsional volume of A^ is zero. We 
can ignore such a sub-domain, since it does not contribute to the integration. 
Applying the above discussion to Ui and J-i for P, wc obtain: 



Gi 



E E 

z' 

[ 



bez"i b'ez^i 
X [l + 



ye{yeAtb')e-^-''+^'+^'''y^ 







E ' 

cez"i\{b} 



-(c-b,y) 



7 r 



aw + 



E 

dez^i\{b'] 



-ic'-b',y) 



(29) 



where 



A, 



66' 



A^' n aJ'. 



(30) 



Our first and second problems are solved when we find how to construct Att' (Vfe G 
V6' G Z^O- 

These problems are interpreted in the words of convex geometry. We intro- 
duce several terms and notations according to the appendix of Ref. [H]. The 
convex polyhedral cone for a finite set S in rt-dimensional Euclidean space is 
defined by: 



C{S) -.^ {V rvv G R" I r„ G R>o, yv G S}. 



(31) 



It is noted that a line is a convex polyhedral cone for S = {v, —v} with a fixed 
vector V. Similarly, it is easy to see that an m-dimensional linear subspace 
R™ (to < n) is also a convex polyhedral cone. 



2 A\_B = {a G A I a B}. Do not confuse with Minkowsky sum A- B = A + {-B) defined 
by Eq. ^ 
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The dual cone of a convex polyhedral cone C is defined by : 

C"" := {ye«:"-\{v,y)>0,\/veC}. (32) 

For a finite set S, C(S')^ is also a convex polyhedral cone. 

Let be the set of points translated from such that point b is moved 
to the origin: 

Zl^ := {c-be Z^-i I c e Z^}. (33) 
Comparing with Eqs. (P7)) and we have 

Af := C(Z,^)^nR^-\ (34) 

and 

Atb' = Ciz'^'f nciz^')"" nR^-\ (35) 

Thus what we have to solve is to find how to construct cones {C{Z^^) and 

C{Zf^,')), dual cones (C(Z^') and C(Zjj,') ), the intersections of dual cones 
(Abf,') and to see whether their volume is zero or not. The algorithms to solve 
these problems are developed in combinatorial geometry, which will be discussed 
in the next section. 

Fig. [T] shows a simple example for the case of two variables: 

P{t) = titl+tltl + tltl+tlt2. (36) 

Fig. [TJa) and (b) show Z^ , Z^, respectively. Vectors a, (3 and 7 satisfies the 
condition (a, v — b) > 0, {(i ,v — b) > 0, {P,v — c) > and (7, i; — c) > for all 
V e Z^. Fig. [UJc) shows the resultant sub-domains. 



Y 




Figure 1: An example of cone and dual cone 



Since these cones are in general still too complex for the integration, we 
divide them into a set of simpler objects. In the case of polytope (bounded 
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polyhedron), it can be decomposed as a union of simplicies, which is called 
triangulation . The algorithms of triangulation is also studied in combinatorial 
geometry. Similarly, a cone can be decomposed as a union of simplicial cones 
as discussed in the next section. 

An (N — l)-dimensional cone C(V) is simplicial when V = {vi, • • • , wjv-i} 
is a set of {N — 1) hnear independent vectors. Let Sbb' be the set of {N — 1) 
vectors characterizing simplicial cones such that 



Abb'= U C{V) 



(37) 



represents a triangulation of Atf ■ With this triangulation, Eq. ([^ is expressed 
by: 



G-E E E 

hgz"; b'ez^i ve5,, 



d'^-^y e{y e C{V)) e-('''+T''+'3''''y) 



ab' 



are 

cez^i\{b'} 



-(c-b',a) 



/3 



(38) 



cez"i\{b} 

A point y = (yi) G C'(y) is parameterized by barycentric coordinate: 

N-l 



Vi 



E(^^)» 



(39) 



Since {vj} is a set of {N — 1) linear independent vectors, the correspondence 
between {yi} and {uj} is one-to-one. We regard V = {vi,--- ,vn-i) being a 
{N — 1) X {N — 1) matrix made arranging column vectors {vj}. Point y is ex- 
pressed hy y — Vu in this matrix notation. We subsequently change the variable 



from ?; to u then from u to z, where Zj 
finally obtain: 

Gi = 



Jacobian is | det / H , ^j- We 



„1 N-l 

E E E idet^i/ d--^.Yi^i 

b<£Z"i b'dz^i ves^b' ° j=i 



(y' +')b+pb' .Vj)-1 



N-l 



1 



E n 



{c-b,Vj) 



(40) 



ab' 



cG2"i\{&} J = l 

The original variables are expressed by: 



E « n 



(c-6',Dj) 



-V] 



n 



(41) 



with Jacobian | det V\Wj ij/Ilfe -^fe • 

As it will be shown in the next section, we can take vectors {vj } on the integer 
lattice Vj e Z^o"^\{0}. Since (c - 6,?;) > for aU y G Af and (c - fo) G Z^^\ 
(c — b^Vj) is non-negative integer. Thus the sub-expressions in the brackets of 
Eq. (HQ]) are polynomials of z. 

For fixed fc, let us consider the scale transformation of vector and variable 



Zk defined by Vk \kVk and Zk 



for Afc > 0. It is easy to see that 
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integration G; is invariant under this transformation. We select Vk on the integer 
lattice as nearest possible to the origin. 

Now all the problems are converted to the following ones in convex and 
combinatorial geometry: 

1. to construct intersections of cones Abi,i — C{Z'^^) r\C(Z^') DW^q^ for 
all b e Z^' and b' e Z^' . 

2. to triangulate Ahb' into Llv£S^^,C{V)- 



4 Algorithms 

4.1 Construction of cones 

First we rewrite A^f,/ into a different form. The following relations hold for 
convex polyhedral cones C, C'l and C2, and for sets of finite points Si and S2 
(see. Ref. [TQ Appendix): 

C{SiUS2) = C{Si) + C{S2), (42) 
C^'' = C, (43) 
C/nCa'' - (C1+C2)'', (44) 

where addition '+' represents Minkowsky sum defined by: 

Ci +C2 :- {a + 6 I a e Ci, 6 e C2}. (45) 

Let be the standard basis of Euclidean space R™, that is 

i;":={ei,e2,--- ,e™}, (46) 
ei = (1,0,0,--- ,0), 62 = (0,1, 0,--- ,0), ••• . (47) 

Set M>o represented by: 

Rt^o i^>ay = C(£;™) C(£;™)''. (48) 
With these relations, Eq. (|55|) is rewritten to: 

= {C{Z'^') + C{Z^') + C{E^-^)f (49) 
= C{Z^' U Z^' U E^-^)"" . 



Thus we obtain 



where 



Aw = CiZbb'V, (50) 



Zbb' ■■= Z'^'UZ^'UE^-K (51) 

What we should know is how to construct a convex polyhedral cone for a given 
finite set and how to construct its dual cone. 
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Let us introduce some terms used in convex geometry. A convex hull of a 
subset T of M™ is the smallest convex set including T and is denoted by conv(T). 
When T is a finite set of half-lines starting from the origin, conv(T) is a convex 
polyhedral cone. When T is a finite set of points, conv(T) is a polytope. When 
T is equal to the set corresponding to a polynomial P, conv(Z^) is called 
Newton polytope of P, which plays important role for analyzing properties of a 
multi-variate polynomial. 

Let Khea, fc-dimensional polytope or convex polyhedral cone in m-dimensional 
Euclidean space. A hyperplane H is called to support a set K at a point x G K 
if x G _ff and K is included in one of the two closed half-spaces limited by 
H. That is, when a vector u is normal to H, either {u,y — x) > 0{yy G K) 
or {u,y — x) < 0(Vy G K) is true. A subset P of if is called a faces of K if 
either P = or P = iiT, or if there exists a supporting hyperplane H oi K 
such that F = H n K. A face of convex polyhedral cone or polytope is also 
a convex polyhedral cone or a polytope, respectively. A face is called a ver- 
tex, an edge or a facet when it is 0-, 1- or {k — l)-dimensional, respectively. 
The set of j-dimensional faces is denoted by Pj (K) and the set of all faces by 
P(Pr) := {0} U*^=o the case of a convex polyhedral cone, Pq is {0} if 

Po ^ 0. It is noted that some special convex polyhedral cones, such as lines or 
the whole space, have no vertices. A face of a polytope is expressed by a set of 
vertices included in the face. On the other hand, an edge of convex polyhedral 
cone is expressed by a non-zero vector laying on the edge. A higher than one- 
dimensional face is expressed by a set of edges included in the face. Let / and g 
be two faces of K. When / is a face of g, this relation is denoted hy f < g. The 
set P(A') forms an abstract complex structure in terms of this binary relation 
-< (see. Ref. [15| Proposition A. 5 and A. 16). 

The convex hull problem of a polytope is to find P(conv(5')) for a set of finite 
points S. This problem is popular and well studied in computational geometry 
(see, for example, Ref. [T^). Several software packages have also been developed 
to solve these problems [T71 [TH] . 

Our problem is to construct an algorithm to find P(C(5)) for a convex 
polyhedral cone C{S). There will be two approaches. One is to modify directly 
an algorithm for a polytope to one for a convex polyhedral cone. It is not so 
difficult and we have prepared one for our test program based on the algorithm 
described in Ref. [IB]. The other is to use the output of algorithms for polytopes. 
We discuss here the latter approach. 

Let us consider a correspondence between a convex polyhedral cone C{S) 
and a polytope conv(S') defined for the same set of finite points S C M™. The 
polytope conv(5') is expressed by: 



Since our set of points Zf,b' includes the origin 0, we assume that G S'. Let 
P(conv(S')) be the set of faces obtained by a convex hull algorithm for a poly- 
tope. Comparing (pij) with ((S^ . C{S) is expressed by: 



conv(S') = A„M I A„ > 0, ^ A„ = 1}. 



(52) 



ties ues 




(53) 



t'EJ^o(conv(5)) 
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Although S includes the origin, it is possible that the origin is not a vertex 
of conv(S'). In this case, there exists a line segment L passing through the 
origin and being included in both conv(S') and C{S). Let the line segment 
L be expressed by L = {Xv | a < A < 6} with some non-zero vector v and 
real numbers a < and 6 > 0. For any y £ (7(5*)^, condition X{v,y) > 
is satisfied for both positive and negative value of A, and then {v, y) = 0. So 
C(S')^ is included in a hyperplane perpendicular to vector v. It means that the 
dimension of C(S')^ is less than m, and this sub-domain does not contribute to 
the integration in our problem. 

Now we assume the origin is a vertex of conv(S'). Then we can assume 
the origin is also a vertex of C{S). If not so, any hyperplane passing through 
the origin does not support C{S) nor conv(5'), which is a contradiction to the 
first assumption. These convex polyhedral cones are called strongly convex. Let 
Tj{C(Sy) be the set of j-dimensional faces of C{S) and J^j(conv(S')) be the 
set of j-dimensional faces of conv(5) including the origin : ^j(conv(5)) = {/ G 
J^j (conv(5)) I e /}. Let a face / e JPj(conv(5')) be defined by / = i?nconv(5) 
with a supporting half-plane H . As H also supports C{S) at the origin, a subset 
g = H n C{S) of C{S) is a j-dimensional face of C{S). It is easy to show that 
this correspondence between a face in Tj{conv{S)) and one in Tj{C{S)) is one- 
to-one and preserves the relation -<. It implies that the structure of abstract 
complex oi J-{C{S)) is the same as one oi J-'{comr{S)). As the result, we obtain 
the following algorithm of constructing C{S) in R™: 

1. Construct J-'(conv(S')) from 5 by a convex hull algorithm for a polytope. 

2. If ^ JFo(conv(5)) then we ignore C{SY , since it does not contribute to 
the integration. 

3. Construct a subset of faces T' {comf{S)) selecting ones including the origin. 

4. Interpret T' {co\w{S)) as the set of faces !F{C{S)). 

The dual cone of C{S) is constructed from T{C{S)) by the dual operation 
shown by Proposition A. 6 of Ref. [T5]. In our case, it is sufficient to know the 
set of edges Ti{C{Sy), since C{S) is expressed by: 

C{SY = { ^iuu\^iue^>o.'iueTl{C{sY)}, (54) 

where u € J-i{C{SY) means to select a non-zero vector u for each edge such 
that u lies on the edge. If dimC(S') = m, an edge of C{SY is obtained as a 
half-line starting from the origin and parallel to the normal vector of a facet of 
C{S). If dimC(5') < m, C(S')^ is not strongly convex and the construction of 
^i(C(S')^) is slightly complicated. However, we do not need to consider this 
discussed before. 

We show that when S" is a set of integer vectors, a vector u G JFi(C(5)^) 
can be selected as an integer vector. Since vertices of conv(S') are included 
in the original set of points S, we can select integer vectors to represent the 
edges of C{S). A facet of C{S) is a subset of hyperplane spanned by some of 
these vectors. A vector representing an edge of C{S)^ is a vector normal to the 
facet. The normal vector is obtained by solving a system of homogeneous linear 
equations. The vector u is solved as a vector w multiplied by the inverse matrix 
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of an integer matrix M . The inverse matrix is constructed by the minor of M 
divided by the determinant det M , where the elements of the minor are integers. 
By adjusting the normahzation of vector ly, we can cancel the denominator 
det M. In this way, we can take v as an integer vector. 

As the result, A;,b' can be expressed with a set of integer vectors {ui , • • • , } C 
Z^-i by: 

Afcb' = I r, > 0}. (55) 

i 

With the fact that our set of points includes E^^^ as a subset, we can 
accelerate the algorithm in the following cases: 

• Case where there exists a non-zero vector v G Z^^ UZ^' and v £ M.^q^. In 
this case, negated vector —v is expressed by a linear combination of vectors 
inE^-^ with non-negative coefficients. It means that C{Z^'UZ^'UE^-^) 
includes line segment connecting v and —v and thus the origin is not the 
vertex. So we can discard such a configuration. 

• Case where there exists a non-zero vector v e Z^' U Z^^ and v e R>o~"'^- 
In this case, vector v is an interior point of C{Z^' U Z^' U E^~^) and is 
not on the edge of the cone. So we can eliminate v from the input set of 
points. That is, the input set of points Z^' U Zj^' U E^^^ can be replaced 
by {(Z^' U Z^')\IK>o~^} U {0} U E^'^. This replacement decreases the 
number of input points, which accelerate the convex hull algorithm. 

These conditions are easily checked before applying the convex hull algorithm. 



4.2 Triangulation of convex polyhedral cones 

We discuss here how to triangulate convex polyhedral cones. Since a convex 
polyhedral cone is not bounded, we convert it to a bounded polytope in order 
to make the problem easier to handle. Let if be a hyperplane perpendicular 
to vector I := (1, 1, ■ • • , 1) and passing through a point wl for w G M>o. As 
C M>(7\ this cone can be truncated to a polytope in cutting by H . Then 
the resulting polytope A[,j,, is expressed by: 

Kb' - {E ^^"^ I < < 1, E < 1}' (56) 

] j 

v'j := Vjw/{I,Vj), {l,Vj)^w. (57) 

Since Vj is an integer vector, (J,Vj) is integer. By selecting w as the least 
common multiplier lcm{(I, Vj) \ Vj}, we can take all of {v^} as integer vectors. 

The triangulation of {N — l)-dimensional polytope is not evident. We re- 
strict ourselves to triangulation methods such that no new vertices are added. 
Even with this restriction, the number of resulting simplicies is not unique. The 
following example (Fig. [5]) shows the case where the different numbers of sim- 
plicies are produced by different triangulations. Let A, B, C, E be points in 
the 3-dimensional space defined by: 

A = (0,0,0), B = (1,0,1), C - (1,1,1), 
D = (0,1,1), E ^ (0,0,2). 
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(a) (b) (c) 

Figure 2: An example of different triangulations 



The convex hull of these points is shown by Fig.[2I^a). In cutting the convex hull 
by a triangle {A, C, E), we obtain two simplicies (Fig.mb)). On the other hand, 
there are three simplicies when it is cut by three triangles {A, B, D), {E, B, D) 
and {B,C,D) (Fig. Hfc)). We do not know triangulation algorithms which 
produce the minimum number of simplicies. We use the following algorithm. 

1. If the number of vertices of the input d-dimensional polytope P is c? + 1 
then P is a simplex. The algorithm returns the set {P}. This condition 
is always satisfied when d = 1 (P is a line segment). 

2. Pick up one vertex u from P. Remove u and all interior points of all 
line segments included in P and starting from u. The resulting geometric 
object is the union of polytopes Qj (j = 1, ■ ■ ■ ,m) of dimension d — 1, 
which are (d — l)-diniensional faces not including u. 

3. For all J = 1, • • • , m, apply triangulation algorithm to Qj recursively, and 
set Sj of simplicies is obtained. 

4. Make the union S = ^jSj. 

5. Construct a d-dimensional simplicies by adding u to each (d— l)-dimensional 
simplex in S. The algorithm returns the set of all constructed simplicies. 

It is noted that this algorithm never falls into the infinite loop. Since step 2 
removes at least one face, the number of recusive applications of the algorithm 
is bounded by the number of faces of the polytope. The number of resulting 
simplicies depends on how vertex u is selected in step 2. One choice is to select 
u when it attaches the maximum number of edges. This condition produces 
relatively smaller number of resulting simplicies, which result will be shown in 
the next section. 

For the triangulation of a cone, one must select the origin for u in the step 2 
at the first time before calling subsequent recursive procedure. This selection 
guarantees that each simplex includes the origin as one of its vertices. Let s„ 
be one of the obtained simplicies and characterized by a set of — 1 linearly 
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independent vectors {vj} as: 

s„ = {J2r,v, I ^r, < 1, < r, (Vj)}. (59) 

i j 

From this simplex, a simplicial cone s'^ is obtained by extending the range of 
the values of {vj}: 

= {T.'^m\^^'-^'y^)}- (60) 

It is noted that vector Vj is an integer vector, for it is one of the vertices in the 
polytope before triangulation. 

5 A test implementation 

We have prepared a program in order to confirm the correctness of our method. 
For the convex hull method we adopted incremental algorithm described in 
Ref. [16] and modified it for convex polyhedral cones. In order to avoid the 
accuracy problem of numerical data, all numerical operation are performed in 
integer arithmetic. We adopted program language python which is equipped 
with arbitrary long integer arithmetic. Input data is prepared by the program 
package described in Ref. |10| . The output integrand of our program is passed 
to the same program package, which continues the subsequent divergence sepa- 
ration and integration steps. 

The procedure of constructing convex hull is compared with qhull package 
[17j . The triangulation procedure is checked by computing the volume of result- 
ing simplicies. Integrated values are compared with ones in published articles 
in several cases. We show the number of decomposed sectors in Table [T] Col- 
umn "H" is cited from Ref. and columns "A", "B", "C", "S" and "X" from 
Ref. [131 



Table 1: The number of sectors 



Diagram 


A 


B 


C 


S 


X 


H 


This 
















method 


Bubble 


2 


2 


2 




2 




2 


Triangle 


3 


3 


3 




3 




3 


Box 


12 


12 


12 


12 


12 




12 


Tbubble 


58 


48 


48 




48 




48 


Double box, — 


775 


586 


586 


362 


293 


282 


266 


Double box, pi 












197 


186 


Double box, Pi — 


1138 


698 


698 




395 




360 


nonplanar 
















D420 


8898 


564 


564 


180 


F 




168 


3 loop vertex (A8) 












684 


684 


Triple box 


M 


114256 


114256 


22657 


10155 




6568 



This test program is too slow for the practical purposes. We are planning 
to prepare a practical implementation in C++. 
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6 Conclusion 



We proposed a new method of IR factorization in sector decomposition em- 
ploying a geometric interpretation of the problem. The original problem is 
converted into a set of problems in convex geometry: (1) construction of in- 
tersection among dual cones of convex polyhedral cones corresponding to each 
polynomial in the original problem and (2) triangulation of them. They are 
solved by the algorithms in combinatorial geometry. This is a deterministic 
method and never falls into a infinite loop. The number of resulting sectors 
depends on the algorithm of triangulation. Our test implementation of this 
method shows the smaller number of sectors compared to other methods based 
on iterated decomposition. 
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